R包limma作差异基因分析

R包limma差异基因分析

limma包最初设计用于分析芯片数据，后经扩展也可适用于RNA-seq数据。通过将read count数据转换为log2-counts per million（logCPM），估计均值-方差（mean-variance）关系以确定在线性建模之前每次观察的权重，之后将数据应用于线性建模，并通过经验贝叶斯统计差异基因。

本篇简介R语言limma包的差异基因分析流程。

下文中使用的示例数据及相关R代码的百度盘链接：

https://pan.baidu.com/s/1WWuJ6aE25QGXpfFzbU91VA

limma差异基因分析

R包安装

limma直接使用Bioconductor安装。

#Bioconductor 安装 limma
BiocManager::install('limma')

library(limma)

准备数据

网盘示例数据“gene.txt”，是一个基因表达数据矩阵，包含10列样本（可分为control和treat两组，每组各5个样本）共计1000行基因。

#基因表达矩阵
exprSet <- read.delim('gene.txt', row.names = 1, sep = '\t', stringsAsFactors = FALSE, check.names = FALSE)
head(exprSet)

limma分析需要指定一试验设计矩阵，即样本分组，要保证基因表达矩阵中的样本名称顺序和分组顺序是一一对应的。

#试验设计矩阵
#注意要保证表达矩阵中的样本顺序和这里的分组顺序是一一对应的
group <- factor(rep(c('control', 'treat'), each = 5), levels = c('control', 'treat'))
design <- model.matrix(~0+group)
colnames(design) = levels(factor(group))
rownames(design) = colnames(group)
design

基因表达值标准化（voom标准化）

计数数据（read count）存在不可忽视的均值-方差关系，原始计数表现出随着计数大小的增加而增加的方差，而对数计数则通常表现为均值-方差趋势减小。

limma包中voom()函数用于标准化RNA-Seq或芯片数据。该函数将计数（所有计数加0.5以避免对数取零）转换为logCPM值，之后将logCPM值矩阵进行标准化。通过估计对数计数的均值-方差趋势，然后根据其预测方差为每个观察值分配权重，用于在线性建模过程中使用权重来调整异方差。

#voom 标准化，详见 ?norm
norm <- voom(exprSet, design, plot = TRUE)

数据标准化后，根据每个基因的平均log2CPM构建线性模型，计算残差，并拟合log2CPM与残差标准差平方根的关系，得到的趋势线（红色平滑曲线）可为样本和基因分配权重。

对于该趋势线，若左侧0起点处所示残差标准差明显偏高，或者0起点处出现上式趋势，则表明数据中存在较多的低表达（low counts数）基因。若觉得它们可能会对后续的计算带来较大的误差，不妨在标准化前手动过滤下。

数据执行标准化后，原本离散程度较大的reads counts矩阵的离散程度降低。

#标准化前后数据比较
par(mfrow = c(2, 2))
boxplot(exprSet)
plotDensities(exprSet)
boxplot(norm$E)
plotDensities(norm$E)

此外，对于voom()标准化功能，除了像上述这样指定原始的基因表达矩阵外，也可将其应用于edgeR的DGEList对象中实现数据标准化，因此limma也常和edgeR结合使用。前文在简介edgeR时，在最后就展示了一例voom()标准化DGEList的例子。

差异分析

接下来就是使用标准化后的基因表达数据，通过加权或广义最小二乘法对每个基因拟合线性模型，并通过经验贝叶斯计算出适度的t统计值、F统计值和差异表达值，最后获得差异基因分析结果。

#线性拟合，详见 ?lmFit
fit <- lmFit(norm, design, method = 'ls')

#确定比较的两组
#后续将计算标记为 1 的组相对于 -1 的组，基因表达值的上调/下调状态
contrast <- makeContrasts('treat-control', levels = design)
contrast

#使用经验贝叶斯模型拟合标准误差，详见 ?eBayes
fit2 <- contrasts.fit(fit, contrast)
fit2 <- eBayes(fit2)

qqt(fit2$t, df = fit2$df.prior+fit2$df.residual, pch = 16, cex = 0.2)
abline(0,1)

#p 值校正、提取差异分析结果，详见 ?topTable
diff_gene <- topTable(fit2, number = Inf, adjust.method = 'fdr')
head(diff_gene, 10)

write.table(diff_gene, 'gene_diff.txt', col.names = NA, sep = '\t', quote = FALSE)

后续自定义根据log2FC值以及校正后的p值等，自定义筛选差异基因就可以了。

火山图示例。

#例如这里根据 |log2FC| >= 1 & FDR p-value < 0.01 定义“差异”
diff_gene[which(diff_gene$logFC >= 1 & diff_gene$adj.P.Val < 0.01),'sig'] <- 'red'
diff_gene[which(diff_gene$logFC <= -1 & diff_gene$adj.P.Val < 0.01),'sig'] <- 'blue'
diff_gene[which(abs(diff_gene$logFC) < 1 | diff_gene$adj.P.Val >= 0.01),'sig'] <- 'gray'

log2FoldChange <- diff_gene$logFC
FDR <- diff_gene$adj.P.Val
plot(log2FoldChange, -log10(FDR), pch = 20, col = diff_gene$sig)
abline(v = 1, lty = 2)
abline(v = -1, lty = 2)
abline(h = -log(0.01, 10), lty = 2)

友情链接